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Abstract 


Contact resistance between the bipolar plate (BPP) and the gas diffusion layer (GDL) in a proton exchange membrane (PEM) fuel cell constitutes 
a significant portion of the overall fuel cell electrical resistance under the normal operation conditions. Most current methods for contact resistance 
estimation are experimental and there is a lack of well developed theoretical methods. A micro-scale numerical model is developed to predict the 
electrical contact resistance between BPP and GDL by simulating the BPP surface topology and GDL structure and numerically determining the 
status for each contact spot. The total resistance and pressure are obtained by considering all contact spots as resistances in parallel and summing 
the results together. This model shows good agreements with experimental results. Influences of BPP surface roughness parameters on contact 
resistance are also studied. This model is beneficial in understanding the contact behavior between BPP and GDL and can be integrated with other 


fuel cell simulations to predict the overall performance of PEM fuel cells. 


© 2006 Published by Elsevier B.V. 
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1. Introduction 


Fuel cells are a promising power technology with a wide 
variety of potential applications. Particularly, proton exchange 
membrane (PEM) fuel cells have received broad attentions due to 
their low operation temperature, low emission and quick startup. 
One of the key technical barriers to the commercialization of 
PEM fuel cells is the cost-effective manufacturing and preci- 
sion assembly of fuel cell stacks to achieve the desired perfor- 
mance. One of such performance measures is the cell potential, 
which decreases from its equilibrium potential during operation 
because of irreversible losses caused by activation, concentra- 
tion and ohmic resistances. Among all these resistances, ohmic 
resistance is dominant under the normal fuel cell operation con- 
ditions. Contact resistance constitutes a significant part of the 
ohmic resistance, especially when stainless steel, titanium or 
molded graphite is chosen as the BPP material [1,2]. 

Contact resistance occurs at all interfaces inside the fuel cell, 
the most important one being at the interface between BPP and 


* Corresponding author. Tel.: +1 734 615 4315; fax: +1 734 647 7303. 
E-mail address: jackhu@umich.edu (S.J. Hu). 


0378-7753/$ — see front matter © 2006 Published by Elsevier B.V. 
doi: 10.1016/j.jpowsour.2006.09.019 


GDL, as shown in Fig. 1. Contact resistance is determined by 
the material properties, surface topology, assembly pressure and 
operation conditions. During assembly, an optimal assembly 
pressure is needed to balance the contact resistance and flow 
resistance in GDL [3]. A high assembly pressure can reduce the 
contact resistance, but the GDL will be over compressed with 
high stress, which results in increased flow resistance. Thus, 
understanding the contact resistance mechanisms between BPP 
and GDL is important in optimizing clamping pressure as well 
as improving the fuel cell performance. 

Some experimental researches have been conducted on the 
contact resistance in PEM fuel cell. Mathias et al. [4] showed that 
contact resistance between GDL and BPP is greater than the bulk 
resistance of GDL or BPP. Ihonen et al. [5] developed a novel 
PEM fuel cell assembly to measure the clamping pressure and 
contact resistances simultaneously for laboratory investigations. 
Results showed that contact resistances depended on clamping 
pressure, gas pressure, current density and temperature. Also, the 
contact resistances of stainless steel could be drastically reduced 
by surface treatments. Lee et al. [6] measured the PEM fuel cell 
performance with a variety of commercially available GDLs 
under various assembly pressures. Each GDL exhibited its own 
optimal assembly pressure due to the differences in mechanical 
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Fig. 1. Schematic structure of a PEM fuel cell. 


properties and porous characteristics, resulting in different con- 
tact resistances in PEM fuel cells. 

All of the above-mentioned studies focused on obtaining the 
contact resistance experimentally. Only a few attempts were 
made on the development of theoretical models of the contact 
resistance in PEM fuel cells. Mishra et al. [7] used a fractal based 
model to predict the contact resistance between GDL and BPP 
and measured the contact resistance experimentally. However, 
the GDL surface roughness parameters, which are important 
inputs for the fractal model, change during compression and 
are difficult to characterize. In a recent work, Zhang et al. [8] 
developed simple computational methods for estimating con- 
tact resistance between BBP and GDL based on experimentally 
obtained constitutive resistance—pressure relations. 

Despite the lack of theoretical models of contact resistance in 
PEM fuel cells, a significant amount of literature exists in mod- 
eling of electrical contact resistance between contacting bodies. 
Most of these models incorporate the contact behavior of a sin- 
gle spherical asperity into a statistical model of multi-asperity 
contact [9-12]. The most recognized one is the Greenwood and 
Williamson (G&W) statistical model [9], which is based on the 
Hertz solution for individual elastic contacts and assumes that 
only asperities originally higher than the separation of the sur- 
faces are in contact. This statistical method accounts for the 
stochastic nature of the interfacial phenomena and has been 
widely used to predict the contact of rough surfaces. 

This paper develops a micro-scale numerical model to predict 
the contact resistance between BPP and GDL in PEM fuel cells 
taking the material properties, surface roughness and clamp- 
ing pressure as inputs. The classical G&W contact model of 
rough surfaces is adapted. Experimental measurements of con- 
tact resistance between BPP and GDL are also conducted to 
validate the model. The remainder of the paper is organized as 
follows: Section 2 introduces the micro-scale contact resistance 
model, Section 3 describes the numerical example and experi- 
ments, Section 4 presents the results and discussions and Section 
5 draws the conclusions. 


2. Micro-scale contact resistance model 


In PEM fuel cells, the BPP surface is rough in nature while 
GDL is a porous medium consisting of randomly distributed 
fibers. Real contact occurs microscopically between BPP asper- 
ities and GDL fibers. Therefore, the topologies of both BPP and 
GDL are important in understanding their interfacial contact 
behaviors. The micro-scale contact model is developed using 


the following procedure: 


(i) The BPP surface topology is simulated as randomly dis- 
tributed asperities based on measured surface roughness 
using profilometrical measurements. 

(ii) The GDL is modeled as randomly distributed cylindrical 
fibers with its total fiber length estimated from the GDL 
porosity and measured fiber diameter. 

(iii) Given a nominal separation between the BPP and GDL, 
BPP asperities in contact with GDL are determined numer- 
ically. 

(iv) The contact area, force and contact resistance of every sin- 
gle contact spot between BPP and GDL are calculated using 
the Hertz theory. 

(v) The total contact resistance is calculated by considering all 
contact spots as resistances in parallel and the total clamp- 
ing force is the summation of the forces on all contact spots. 

(vi) Experimental measurements of contact resistance were 
conducted to validate the modeling results. 


2.1. BPP surface topology simulation 


All surfaces are inherently rough. The surface of a BPP, 
whether molded using graphite or formed with stainless steel, 
contains surface roughness, which determines the contact behav- 
ior. Consistent with the classical statistical contact models of 
rough surfaces, the BPP surface is assumed to be covered with 
asperities whose summits are all spherical in shape with the same 
radius R;. The summit height follows a normal distribution. The 
summits are also assumed to be uniformly distributed spatially 
with a known density Dsum, measured in “number of summits per 
unit area”. Three parameters, summit radius Rj, standard devi- 
ation of summit height os and summit density Dsum are needed 
to describe the surface roughness. According to McCool [12], a 
non-dimensional parameter œ is introduced, 


eo (1) 


where ø, o’. and o” are the root mean square of surface height, 
slope and second derivative of a surface profile, respectively. 
The mean summit radius is expressed as 


3771/2 


I 
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and the variance of the summit height distribution can be calcu- 
lated from o [11], 


pe (i u -2 B 


Q 


In a two-dimensional surface profile, the local highest point 
is a peak. The peak density is 


Dpeak =~ (4) 


The surface profile was obtained using a profilometer with 
a lateral resolution of 0.5 um. Several scans in different parts 
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Fig. 2. The generated BPP surface (unit: mm). 


of the BPP surface were conduced to obtain one set of rough- 
ness parameters for the entire surface. The surface roughness 
parameters obtained from the average values of several scans 
are Dpeak = 98 mm~!, Ry = 3.67 um and ø; =3.55 pm. 

A three-dimensional surface profile has identical statistical 
characteristics in every two-dimensional direction. Therefore, 
the surface summit density is assumed to be Dies although it 


was shown to be slightly larger than Dik [12]: 

According to these surface roughness parameters, a surface 
profile is generated to simulate the BPP rough surface. A surface 
is generated with 98 x 98/mm? randomly distributed spherical 
summits with 3.67 um in radius. The summit height is normally 
distributed with a standard deviation of 3.55 wm. A sample area 
of the generated surface is illustrated in Fig. 2, which is obtained 
from one simulation. 

This generated surface provides the same surface roughness 
characteristics as the measured roughness of the real BPP. It also 
gives the position and height of each asperity summit, which are 
important inputs for the numerical contact resistance model. 


2.2. GDL structure simulation 


The GDL is made of carbon fiber paper or carbon fiber cloth. 
The carbon fiber paper is one of the primary materials due to 
its high porosity (>70%) and good electrical conductivity. It is 
made from polyacrylonitrile-precursed-carbon fiber, the same 
material as used for reinforced composite. During the manufac- 
turing process, the chopped carbon fibers are dispersed in water 
with binders and dried layer by layer to achieve the required 
thickness. Carbon fibers with a diameter of approximately 7 ym 


and different lengths are randomly distributed to form the carbon 
fiber paper [4], which can be seen in Fig. 3. 

As shown in the SEM micrographs, the GDL used in this 
study can be approximately treated as a layered structure with 
binders between fibers. The binder thickness between layers 
varies at different locations of GDL, from close zero to as much 
as 6 um. The majority of the binder thickness is about 4 wm. 
According to the images and assumptions, the carbon fiber paper 
can be characterized as follows: 


e the carbon paper is made of multi-layers of carbon fibers; 

e the carbon fiber is cylindrical in shape with a diameter dfbver 
of approximately 7 um; 

e the carbon fibers are randomly distributed in length and ori- 
entation at each layer; 

e the binder thickness ôbinder between two adjacent layers is 
approximately 4 wm. 


The total fiber length in a unit area of this sample can be 
obtained as 


Vap x (1 — €) 


e= 
er — 2 
1/47dfher 


(5) 


where Lif is the total fiber length in the unit area, Vapi the 


volume of the GDL sample and ¢ is the porosity. The fiber length 
in each layer is the total fiber length divided by the number of 
layers. 
tot 
liter = <b © 
ÔGDL/(dfiber + Sbinder) 

Based on these characteristics of the carbon fiber paper, one 
layer of GDL is simulated with randomly distributed carbon 
fibers, as shown in Fig. 4. The GDL structure is generated as 
follows: the location of center point and orientation of each fiber 
varies independently and uniformly in this area. The length of 
each fiber is assumed to be uniformly distributed from 0 to the 
diagonal length of this area. When a fiber intersects with existing 
fibers and boundaries, it is cut at the point of intersection, and 
the remainder of the fiber turns out to be the new fiber length. 
Hence, fibers, which appear late have more chance to be cut and 
become shorter. 

The fiber locations and lengths are deterministic for each 
simulation. Every individual contact spot can be located based 


50 um 


Fig. 3. SEM micrographs of carbon fiber paper. 
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Fig. 4. Simulation of one layer of carbon fiber in a GDL. 


on the information of the relative positions between fibers and 
BPP asperities. 


2.3. Contact resistance numerical model 


The contact resistance between BPP and GDL is governed by 
the surface topography and material properties of the contacting 
pairs. The BPP surface is a rough surface with spherical asper- 
ities, which are in the same order of magnitude as the carbon 
fiber diameter (~7 um). For asperities with heights between 0 
and 20,, the contact with GDL is in the first carbon fiber layer. 
The contact problem is then simplified as asperities contacting 
with one layer of carbon fibers while neglecting carbon fiber 
surface roughness. 

The behavior of an individual point of contact is known from 
Hertzian equations [13]. When a cylinder contacts a sphere 
with nearly the same radius, as in this study, the contact spot 
is close to a circle and the relation between contact area a 
and the load F can be expressed approximately in terms of 
deformation ô as, 


a ~ Rd (7) 


4 
Fx z E” Re" 282 (8) 


Carbon fiber ; . 
Spherical asperity 


and 


—1 
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where E1, E>, v1, v2 are the Young’s moduli and Possion’s ratios 
of the two contacting bodies, respectively. Re is the equivalent 
radius of the principal radii of curvature of the surfaces at the 
contact origin. 


R 
Re Ri (10) 
Ri + R2 


where R2 denotes the carbon fiber radius. 
According to Holm [14], the electrical constriction resistance 
of this single contact was: 


_ +p 
4r 


where r is the radius of the contact area. oq and p2 are the 
resistivities of the two contacting bodies, respectively. 

Three basic assumptions for the contact model are made: (1) 
asperities are far apart and there is no interaction among these 
asperities; (2) there is no bulk deformation in the bipolar plate; 
(3) contact is entirely elastic. 

Fig. 5 shows a two-dimensional illustration of the relative 
distance between asperities and fibers. In reality, the fibers are 
not necessarily parallel and can be located in any direction. O 
and O’ denote the center of spherical asperity and carbon fiber. If 
the distance OO’ is less than the summation of two radii R; + Ro, 
the cylinder and the asperity are in contact. Otherwise, there is no 
contact between them. For a given separation d, the deformation 
between the asperities and fibers is easy to calculate, so the 
contact area, force and resistance for each contact spot can be 
determined. The total contact area is the summation of all contact 
spots and the total contact resistance is calculated by considering 
the resistance of all contact spots in parallel. The total force is 
the summation of all contact forces. 


R (1) 


3. Numerical example and experiments 
3.1. Numerical example 
The above procedure is implemented numerically. A rough 


BPP surface of 4mm x 4mm is simulated at an initial sepa- 
ration of 7.5 um from a GDL layer with the same area. In 


Center line of carbon fibers 


o 
OO'SR +R 


o 
OO`=R+R2 


oO 
OO'>R }+R > 


Lower envelop surface of GDL 
— Summit mean plane of BPP 


Fig. 5. Two-dimensional illustration of relative position between BPP asperities and carbon fibers. 
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Table 1 

Material properties for the BPP and the GDL 

Properties BPP GDL 
Thickness (mm) 5 0.11 

Area (mm?) 101.6 x 101.6 100 x 100 
Porosity 80% 
Young’s modulus (GPa) 10 

Electrical resistivity-through plane (uQ m) 190 800 


Table 2 
Inputs parameters for the numerical contact model 
Parameters Value 
BPP 
Asperity peak density Dpeak (#mm7!) 98 
Summit radius R (um) 3.67 
Summit standard deviation os (jm) 3.55 
GDL 
Fiber diameter dgher (wm) 7 
Total fiber length Li! „ (mm mm”) 572 
Fiber length in one layer lfiber (mm mm~?) 57 


this study, the BPP is a grade FU 4369 graphite plate from 
PEM Technology Inc and the GDL is Toray TGP-H-030 from 
Toray Industries, Inc. All of the relevant material properties are 
listed in Table 1. Table 2 lists the inputs parameters for the 
numerical models. Parameters for BPP are obtained based on 
the profilometerical scans. Parameters for GDL are estimated 
from the GDL porosity and SEM micrographs. Other inputs to 
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the numerical model include carbon fiber material properties. 
The carbon fiber has more favorable mechanical and electrical 
properties in the longitudinal direction than in the transverse 
direction [4]. However, the transverse material properties have a 
more significant influence on the contact resistance in PEM fuel 
cells. Hence, such material properties are used in the numerical 
model. The transverse compressive modulus of Toray carbon 
fiber is 3.2 GPa [15]. The transverse electrical resistivity is 
70 uQ m, which is estimated using the Bruggman correlation 
[16], 


pı = pop. — £)!5 (12) 


where papu is the through plane resistivity of the GDL. 

Based on the simulated BPP surface and GDL structure, the 
relation between clamping pressure and contact resistances can 
be obtained by changing the separation. 


3.2. Experimental validation 


Experimental investigations were conducted to validate the 
numerical model results. Two experimental setups were built 
to measure the contact resistance [4,7,8]. Setup | as shown in 
Fig. 6(a) was built using a stack of a GDL and two graphite 
BPPs. This stack was inserted between two copper plate cur- 
rent collectors. Plexiglas plates were used for insulation. The 
measured resistance from Setup 1 includes the bulk resistances 
of two BPPs, the bulk resistance of GDL, contact resistances 
between cooper plates and BPPs, BPPs and GDL. Setup 2 uses 
a similar stack but with only one BPP between two copper plates 
in order to extract the contact resistance between BPP and GDL. 


Insulation 


Bipolar plate 


Copper plates 


(b) 


Fig. 6. (a) Schematic of two experimental setups and (b) picture of Setup 2. 
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An MTS machine was used to provide the clamping load and a 
DC milli-ohmmeter (GW-Instek GOM-802) was used to mea- 
sure the resistance with a resolution of 0.1 pQ. 
The contact resistance between BPP and GDL can be deduced 
as [7] 
Res; — Res2 — Rapp — Rept 


Reontact = 7 ( 14) 


where Res, and Res» are measured resistances from Setups 1 and 
2, respectively. Rppp is the bulk resistance of graphite BPP and 
Rept is the bulk resistance of GDL. Rgpp and Rapu are calcu- 
lated according to their bulk resistivities. The change of the bulk 
resistance of BPP and GDL during compression is neglected. A 
series of compression pressures from 0.5 to 3 MPa were applied 
and the corresponding contact resistances were measured. Under 
each clamping pressure, the contact resistance measurements 
were repeated four times to obtain the average values. Two GDL 
samples are used and results are identified as Experiments 1 
and 2. 


4. Results and discussions 


Results from the numerical model and experiments are pre- 
sented in Fig. 7. Results from Experiments | and 2 are very 
comparable. The discrepancy between two experiments is less 
than 5%, which illustrates that the contact resistance measure- 
ment is repeatable. For every clamping pressure, the simulation 
was repeated five times. Results from repeated simulations show 
a small range of variability, in particular, when the clamp- 
ing pressure becomes large. This is because, at high clamping 
pressure, the number of contact spots increases and the calcu- 
lation is more accurate. The maximum relative error among 
different simulation runs is less than 3.5%. Furthermore, the 
numerical prediction shows the same trend as the experimental 
results and the difference is less than 20%. This consistency 
indicates that the numerical model captures reasonably well 
the contact phenomena between BPP and GDL in PEM fuel 
cells. 

The number of contact spots increases rapidly when two con- 
tacting parts start to come into contact. The real contact area 
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Fig. 7. Comparison of experimental data with numerical prediction. 


increases fast at this early stage because many asperities are 
coming into contact with the fibers. Thus, the contact resistance 
changes greatly during the initial clamping pressure increase as 
shown in Fig. 7. As the distance between two surfaces decreases 
further, only a few new asperities are becoming involved in 
the contact. Contact resistance decreases mainly due to the 
area increase of the existing contact spots. This slight increase 
of the contact area results in only small changes in the con- 
tact resistance. Although greater clamping pressure can reduce 
the contact resistance, high pressure may damage BPP, GDL 
and obstruct gas flow. From fuel cell performance prospective, 
electrical resistance and flow resistance need to be optimized 
simultaneously to obtain a proper clamping pressure. 

In the numerical model, several assumptions were made to 
model BPP surface and GDL structure. In modeling the BPP sur- 
face roughness, the assumption of spherical asperities with iden- 
tical radius is consistent with the G&W model. Some researchers 
have developed contact models with different asperity shapes 
and radii, but only shown that the G&W model is nevertheless 
quite accurate [10-12,17]. 

Although the length of each fiber is initially assumed to fol- 
low a uniform distribution, the final distribution will change 
after cutting. However, the different fiber length distributions 
after cutting will not change the total contact resistance. This is 
because BPP surface is a random surface, and every location of 
the GDL surface is statistically equivalent. To verify this numer- 
ically, a simplified GDL structure with the same total fiber length 
is modeled, in which all the fibers are in the horizontal direction. 
The difference between these two results is within 3%, therefore 
this assumption is reasonable. 

The influences of BPP surface roughness parameters are also 
investigated. As shown in Fig. 8, when the peak density increases 
from 50 to 150mm !, contact resistance decreases because 
more asperities are in contact. However, the contact resistance 
decreases about 14%, not as much as the change of peak den- 
sity. For a given distance between two contacting surfaces, the 
larger summit density will result in an increased contact area 
and a smaller contact resistance, but larger clamping pressure is 
needed accordingly. In order to compare the resistances at dif- 
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Fig. 8. Influence of BPP summit density. 
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* R: =3.67um 
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Fig. 9. Influence of BPP summit radius. 
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Fig. 10. Influence of BPP summit standard deviation. 


ferent peak densities with a given pressure, the distance between 
the contact surfaces will be different. That is, smaller distance 
for low density and bigger distance for high density are required, 
respectively. The combined effect of the distance and the peak 
density makes the contact resistance not sensitive to the sum- 
mit density changes. The contact resistance is also insensitive to 
summit radius change, as illustrated in Fig. 9. 

The summit standard deviation of a BPP surface is another 
important surface roughness parameter. As shown in Fig. 10, 
when o; changes from 3.55 to 7 um, the contact resistance 
increases by about 20%. For the surface with os =7 um, the 
likelihood to have high asperities is higher than that of smaller 
standard deviation. This results in the fact that the contact spot 
concentrates in a few asperities at the same clamping pressure. 
Thus, the contact area is smaller and the contact resistance turns 
out to be larger. 

The contact resistance is more sensitive to summit standard 
deviation than summit density and summit radius. Hence, BPP 
surface should be fabricated with small standard deviation to 
ensure the consistency of contact resistance and the performance 
of PEM fuel cells. 


5. Conclusions 


A micro-scale contact model was developed to predict the 
contact resistance between BPP and GDL in PEM fuel cells. 
BPP surface roughness was simulated by adopting the classical 
statistical contact model and the GDL was modeled as randomly 
distributed fibers with estimated total fiber length. According to 
these two simulated contacting surfaces, contact spots between 
BPP and GDL can be determined numerically given a separation 
of these two surfaces. The contact status for every single contact 
spot was calculated using the Hertz theory. The total resis- 
tance and pressure were obtained by summarizing the results 
from each contact spot. Compared with experimental results, 
the modeling results showed good agreements with less than 
20% discrepancy. Influences of BPP surface roughness parame- 
ters on contact resistance were also studied. It was found that the 
summit standard deviation has greater impact than other surface 
roughness parameters. The model developed in this study may 
be applied to predict the contact resistance of GDL in contact 
with other BPP materials, only if the GDL characteristics and 
BPP surface topology are determined. This micro-scale contact 
model is beneficial to understand the basic mechanisms of con- 
tact behavior between the rough surface and a fibrous medium 
and can be integrated with other fuel cell simulations to predict 
the overall fuel cell performance. 
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